package twoD.hofem;

/**
 * Class provides functions for integration using Gaussian quadrature rule.
 * 
 * @see FunctionRToR
 * @see GaussPoints
 * 
 * @author <font style="background-color: black; color: white;">&nbsp;:: Team 1
 *         ::&nbsp;</font><a href="mailto:luowei.de@googlemail.com"><font
 *         style="background-color: red; color: white;">&nbsp; Luo
 *         &nbsp;</font></a><a
 *         href="mailto:enes.siljak@ruhr-uni-bochum.de"><font
 *         style="background-color: green; color: white;">&nbsp; Siljak
 *         &nbsp;</font></a><a href="mailto:engr_asifkhan@yahoo.com"><font
 *         style="background-color: blue; color: white;">&nbsp; Khan
 *         &nbsp;</font></a>
 * 
 */
public class Integrator {

	private static double integratePart(FunctionRToR f, Interval interval, int n) {
		double result = 0, x, ksi, weight;

		double a = interval.getLower();
		double b = interval.getUpper();

		for (int i = 0; i < n; i++) {
			ksi = GaussPoints.getPoint(n, i);
			weight = GaussPoints.getWeight(n, i);

			x = ((b - a) * ksi + b + a) / 2;
			result += weight * f.valueAt(x);
		}

		return result * (b - a) / 2;
	}

	/**
	 * Integrates polynomial with the given number of Gauss points.
	 * 
	 * @param f
	 *            Function to be integrated
	 * @param numbOfGaussPoints
	 *            Number of Gauss points
	 * @return Value of integral of the function <b>f</b> over its domain.
	 */
	public static double integratePolynomial(FunctionRToR f,
			int numbOfGaussPoints) {
		double[] bps = f.getBreakpoints();

		if ( (bps == null) || (bps.length == 0) ) {
			return integratePart(f, f.getDomain(), numbOfGaussPoints);

		} else {
			Interval interval;
			double result;

			// first part
			interval = new Interval(f.getDomain().getLower(), bps[0]);
			result = integratePart(f, interval, numbOfGaussPoints);

			// parts 2 to n-1 defined by break points
			for (int i = 0; i < bps.length - 1; i++) {
				interval = new Interval(bps[i], bps[i + 1]);
				result += integratePart(f, interval, numbOfGaussPoints);
			}

			// last part
			interval = new Interval(bps[bps.length - 1], f.getDomain().getUpper());
			result += integratePart(f, interval, numbOfGaussPoints);

			return result;
		}
	}

	private static double adaptiveIntegration(FunctionRToR f, Interval interval) {

		int numbOfGaussPoints = 10;
		
		Interval leftHalf = new Interval(interval.getLower(), (interval.getLower() + interval.getUpper()) / 2);
		double L = integratePart(f, leftHalf, numbOfGaussPoints);

		Interval rightHalf = new Interval((interval.getLower() + interval.getUpper()) / 2, interval.getUpper());
		double R = integratePart(f, rightHalf, numbOfGaussPoints);

		double LR = integratePart(f, interval, numbOfGaussPoints);
		
//		if (Math.abs(L + R) < 1E-12)
//			return 0;

		double relError = Math.abs((L + R - LR) / (L + R));

		double dx = interval.getUpper() - interval.getLower();
		if ( (relError < 1E-8)  || (dx < 0.01) )
			return L + R;

		else {
			L = adaptiveIntegration(f, leftHalf);
			R = adaptiveIntegration(f, rightHalf);
			return L + R;

		}

	}

	private static double integrateNonpolynomial(FunctionRToR f) {
		double[] bps = f.getBreakpoints();

		if ( (bps == null) || (bps.length == 0) ) {
			return adaptiveIntegration(f, f.getDomain());

		} else {
			Interval interval;
			double result;

			// first part
			interval = new Interval(f.getDomain().getLower(), bps[0]);
			result = adaptiveIntegration(f, interval);

			// parts 2 to n-1 defined by break points
			for (int i = 0; i < bps.length - 1; i++) {
				interval = new Interval(bps[i], bps[i + 1]);
				result += adaptiveIntegration(f, interval);
			}

			// last part
			interval = new Interval(bps[bps.length - 1], f.getDomain()
					.getUpper());
			result += adaptiveIntegration(f, interval);

			return result;
		}

	}

	/**
	 * Integrates function over domain of its definition.
	 * 
	 * <br>
	 * <br>
	 * In case that function is:
	 * <ul>
	 * <li>Polynomial: Function calculates number of integration points as:
	 * <b>Ng = (p + 1) / 2</b>, where <b>p</b> is order of polynomial, and uses
	 * Gaussian quadrature rule.</li>
	 * <li>Non-polynomial: function performs adaptive integration combined with
	 * Gaussian quadrature rule.</li>
	 * </ul>
	 * 
	 * @param f
	 *            Function to be integrated
	 * @return Value of integral of the function <b>f</b> over its domain.
	 * 
	 * 
	 * @author <font style="background-color: black; color: white;">&nbsp;::
	 *         Team 1 ::&nbsp;</font><a
	 *         href="mailto:luowei.de@googlemail.com"><font
	 *         style="background-color: red; color: white;">&nbsp; Luo
	 *         &nbsp;</font></a><a
	 *         href="mailto:enes.siljak@ruhr-uni-bochum.de"><font
	 *         style="background-color: green; color: white;">&nbsp; Siljak
	 *         &nbsp;</font></a><a href="mailto:engr_asifkhan@yahoo.com"><font
	 *         style="background-color: blue; color: white;">&nbsp; Khan
	 *         &nbsp;</font></a>
	 * 
	 * @see GaussPoints
	 * @see FunctionRToR
	 * 
	 * 
	 */
	public static double integrate(FunctionRToR f) {

		if (f.getP() == -1) {
			return integrateNonpolynomial(f);

		} else {
			// p.max.exact = 2*ng - 1
			// ng = (p + 1)/2

			int numbOfGaussPoints = (int) Math.ceil((f.getP() + 1.0) / 2.0);
			if (numbOfGaussPoints > GaussPoints.MAX_N)
				return integrateNonpolynomial(f);
			
			else
				return integratePolynomial(f, numbOfGaussPoints);
		}
	}

	public static double integrate(int np, FunctionR2ToR u) {
		double sum = 0;

		for (int i = 0; i < np; i++) {
			double xi1 = GaussPoints.getPoint(np, i);
			double w1 = GaussPoints.getWeight(np, i);

			for (int j = 0; j < np; j++) {
				double xi2 = GaussPoints.getPoint(np, j);
				double w2 = GaussPoints.getWeight(np, j);

				sum += w1 * w2 * u.valueAt(xi1, xi2);
			}
		}
		return sum;
	}

}
